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ABSTRACT 

Several population synthesis models now predict integrated colors of simple stellar 
populations in the mid-infrared bands. To date, the models have not been extensively 
tested in this wavelength range. In a comparison of the predictions of several recent 
population synthesis models, the integrated colors are found to cover approximately the 
same range but to disagree in detail, for example on the effects of metallicity. To test 
against observational data, globular clusters are used as the closest objects to idealized 
groups of stars with a single age and single metallicity. Using recent mass estimates, we 
have compiled a sample of massive, old globular clusters in M31 which contain enough 
stars to guard against the stochastic effects of small-number statistics, and measured 
their integrated colors in the Spitzer /IRAC bands. Comparison of the cluster photom- 
etry in the IRAC bands with the model predictions shows that the models reproduce 
the cluster colors reasonably well, except for a small (not statistically significant) offset 
in [4.5] — [5.8]. In this color, models without circumstellar dust emission predict bluer 
values than are observed. Model predictions of colors formed from the V band and the 
IRAC 3.6 and 4.5 /im bands are redder than the observed data at high metallicities and 
we discuss several possible explanations. In agreement with model predictions, V — [3.6] 
and V — [4.5] colors are found to have metallicity sensitivity similar to or slightly better 
than V — K s . 

Subject headings: Galaxies: individual: Messier number: M31, Galaxies: star clusters: 
general, Infrared: stars, Stars: Population II 

1. Introduction 

Star clusters are popular probes of the history of galaxy and star formation. As compact sys- 
tems with relatively homogeneous properties, their stellar and dynamical histories can be modeled 
with reasonable accuracy, at least compared to their parent galaxies. Studies of extragalactic star 
clusters have experienced a major step forward with the availability of the Hubble Space Telescope, 
which provided the spatial resolution necessary to identify clusters in distant galaxies and resolve 
those in nearby galaxies into individual stars, and large ground-based telescopes, which provide 
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the light-gathering power to study these faint objects spectroscopically. The study of extragalactic 
GCs has provided important cl ues to the process of ga laxy assembly and complements the direct 
study of high-redshift galaxies (jBrodie fc Straderll2006l ). 



Star clusters also provide important tests of simple stellar population (SSP) models constructed 
via population synthesis methods. Although the evidence is now overwhelming that globular clus- 
ters are not si ngle-age, single-met allicity populations as previously believed (see the review by 
Bragagliall2010l ). they are still less complex than the population mix found in even the smallest galax- 
ies. As such, they provide important obser vational tests for population synthesis models, and there 



have been many examples of such tests ( Peacoc k et al 



2010| ; IPessev et al.1 12008 : iBarmbv fc Huchra 



2000 



2011 



Riffe l _et all l20ld : IConrov k Gunn 



Bica fc Alloinl I1986T). Although becoming in 



creasingly important in the study of nearby stellar populations (e.g. iSheth et al.ll2010i ). and having 
the advantage of reduced extinction sensitivity, the model predictions at mid-infrared wavelengths 
A > 3 /im, have not been well-tested. Mid-infrared observations from ground-based telescopes are 
impractical for all but the brightest sources, due t o the high background, and the existing stel- 



lar libraries may not yet be adequate (although see IVandenbussche et al.1 12002 : ISloan et al 



2003 



Hodge et al.ll2004l ). The ATLAS9 model atmospheres, on which many population synthesis models 
are at least partly based, have line opacity spectra calculated to a maximum wavel ength of 10 /xm. 



Mod eling infrared stellar spectra with a Rayleigh-J eans law is quite co mmon (e.g. IZasowski et al 



2009), but this is "clearly a crude approximation" (jMarigo et al.ll2008l ). Those authors state that 
bolometric corrections for the coolest stars are uncertain by up to a few tenth s of a magnitude. 
In the ir investigation on the use of model atmospheres for infrared calibration, iDecin &: Eriksson 
(|2007l ) found that model spectra are accurate to 3% in the mid-IR but also noted that "for the 
purpose of calculating theoretical spectra in the mid to far-IR, some line lists are still far from 
complete or accurate". 

The Spitzer Space Telescope probes these mid-infrared wavelen gths. While young star c lusters 



have been the subject of numerous investigations with Spitzer (e.g. iGutermuth et al.ll2009l ). there 
has been relatively little analysis of the mid-infrared properties of old star clusters. Old star clusters' 
declining spectral energy distribution in the mid-infrared, in combination with the presumed lack 
of interesting spectral features, is presumably the reason for this lack of attention. Mid-infrared 



observation s have been used in stud i es of mass-losing stars in individual clusters (jBoyer et al 



200fil . bood, bold; lOriglia et all 120071. boioh as well a s searches for cold dust in the intracluster 
medium ( Matsunaga et al. bood : Barmbv et al. 2009 ). One of the few studies of extragalactic 
globular cluster systems is that of ISpitler et al.l (|2008l ). The present paper provides mid- infrared 
photometry of massive globular clusters in the Andromeda galaxy, M31, taken with the Spitzer 
Space Telescope, and compares th em to population sy nthesis models which pred ict magnitudes in 
the IRAC bands, n amely GALEV dKotulla et alJ boogh. FSPS (jConrov et alJbood l and models from 
the Padova group (jMarigo et al.ll2008l ). Only the four IRAC bands are used here . M31 clusters are 
generally too faint to be detected in the MIPS observations of lGordon et al.l ()2005l ). and photometry 
of the handful of Milky Way clusters observed with the MIPS instrument is published elsewhere 
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(|Bover et aU 12000 . 12008 : barmbv et all 12003 ). 

Thr oughout this work we use the Vega magnitude system and assume a distance to M31 of 
780 kpc (jMcConnachie et al.ll2005l ). for which 1" corresponds to 3.78 pc. 



2. Data and Models 



2.1. Cluster Mass Limit 



Star clusters contain a limi t ed number of stars and photometry of star clusters must consider 



stochastic effects (Buzzoni 



1989 



Renzinilll998l ). The integrated light of a cluster can be dominated 



by a few bright s tars, particularly in the in frared where the luminosity of the asymptotic giant 
branch is highest (jFouesneau &: Lanconll2010l ). For this first comparison of star cluster photometry 
with SSP models in the mid-infrared, we wanted to avoid stochastic effects as much as possible, 
so we sought to select only the most massive cl usters. One criterion f or cho osing a selection limit 
is the lowest luminosity limit (LLL) defined by ICerviho &: Luridianal (120041 ). This limit is defined 
for a specific SSP model as the luminosity of the brightest star included in t he model, at the age, 
metallicity, and wavelength under consideration. For 10 Gyr old clusters, ICerviho fc Luridiana 



(|2004j) found the total mass corresponding to the LLL in the K-band to be approximately 10 5 M ; 
those authors showed that a total mass of about 10 times this limit is required in order for stochastic 
effects to cause less than a 10% dispersion in the total luminosity. 

Related to the LLL approach is the concept of the "effective number of stars" emitting at 



a given wavelength, introduced by iBuzzonil (119891 ) . With N e g defined via 5L to t/L 



J tot 



iV 



-1/2 



ell' 



the value of N e q for a 10% uncertainty in luminosity is N e s = (0.1) 



100. For a 10% color 



uncertainty, assuming the same N e g for both bands, N e g > 200 (IBuzzonil 120051 ). The Evolutionary 
Population Synthesis (EPS) website ^ tabulates model predictions of N e g(X) as a function of age for 
populations with L^ Q \ = 1O 1O L . For ages of 10-15 Gyr, typical values of iV e g(4//m) range from 
10 6.6 _ 10 6.85 depending on the metallicity, IMF and horizontal branch treatment; for ages of 5 
Gyr, N e g(AfJLm) ~ 10 6 ' 4 . To scale N e g from Lbol = 1O 1O L to a fixed mass, we used the M/L^ \ 
tabulated by the EPS models. For [M/H] = —1.27, a typical value for the globular clusters under 
consideration, N e g values for a 10 6 M population are 300 (15 Gyr), 270 (10 Gyr), and 260 (8 
Gyr). These suggest that clusters with masses of 10 6 M should not have their integrated mid-IR 
photometry strongly affected by stochastic effects, while clusters with mass 10 5 M would have 
iV e ff well below the requirements. We used this result to set our lower mass limit at 10 6 M . 



^ttp : //www .bo . astro . it/$\sim$eps/home .html 
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2.2. Cluster Sample Selection 



For this first comparison between SSP models and star cluster colors, we wanted to remove 
age as a variable if possible, and select only old clusters. The clusters we considered all belong 
to the Milky Way, M31, or the LMC. While many other Local G roup dwarf galaxies have a few 
globular clusters each (see the compilation by iForbes et al.l 12000) , the few available estimates of 
those clusters' masses in the literature show the largest of them t o have masses about 10 5 M ? , 



M33 has numerous star clusters, with a recent catalogue given by ISaraiedini h Mancond (|2007l ); 
however it has few truly old clusters, and few clusters with mass estimates comparable to those 
available in ot her gal a xies. For Milky Way globular clusters we take the inclusion of a cluster in 
the catalog of iHarria (|2010i ) to mean that its age is > 10 Gyr. For L MC clusters we considered 
only clusters classified as having ages > 10 Gyr bv lPessev et al.l (|2008l ); these ages were compiled 
from the results of color-magnitude diagra m (CMP) fitting. A ge determination for M31 clusters 
is more difficult, and, with one e xception (IBrqwn et al.1 12004! ) . CMD-fitting ages from the main 
sequence turnoff are not available. iPerina et al.l (|201ll ) noted that deriving M31 cluster ages from 
integrated spectroscopy can be problematic. They found that two clusters previously suspected to 
be of intermediate age based on integrated spectra had blue horizontal branches and were likely of 
comparable age to the oldest Galactic GCs. They also found significant differences between ages 
determined from spectral energy distribution (SED) fitting and the CMD analysis. With the caveat 
that ages from spectroscopy may be unreliable, a uniform determination of ages for the M31 GC s 
reduces scatter due to mixing of different methods, and we use the values bv lCaldwell et al.l (|201ll ) 
as the basis for selecting our sample. These were der ived by comparing Lick indices from integrated 
spectroscopy to the SSP models of ISchiavonl (|2007l ). An important limitation of this analysis is 
that ages could not be determined for clusters with metallicities [Fe/H] < — 1. It is possible that 
some younger clusters could be lurking in the sample, although we considered only M31 clusters 
that had not been identified as young by any previous work. 

There are a number of methods for estimating the total mass of a star cluster. While velocity 
dispersions either from integrated spectra or spectra of individual stars provide the most direct mea- 
surements of the gravitational potential of a cluster, such measurements are not available for many 
Local Group star clusters, particularly in M31 which has the largest population. A more widely- 
applicable method involves measuring the total luminosity and multiplying by a mass-to-light ratio. 
Total luminosity can be measured by either large-aperture photometry or fitting and extrapolation 
of a surface br ightness profile . For M 31 clusters, a comprehensive set of mass estimates is available 
in the work of lCaldwell et al.l (120111). who est imated masses by multiplying reddening-corrected V- 
band luminosities (from ICaldwell et aDhood ) by M/Ly = 2. With the goal of selecting a uniform 
sample, we used these mass estimates exclusively to s elect the final M31 sa mple; a few additional 
outlying clusters for which masses were not given by ICaldwell et al.l (120111 ) but which have mass 
estimates in their online catalog (G001, B358, MGC-1, B468, B293) were also i ncluded . For M ilky 
Way and LMC clusters, we used near-infrared i^-band magnitudes as given by ICohenl (|2006l ) and 
Pessev et al.l (|2006l . 120081 ) . The near-infrared is often used as an indicator of stellar mass in studies 
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of external galaxies since it is ex pected to be dominat ed by low-mass stars and is not strongly 
sensitive to recent star formation (IBe 
to luminosities using M &j ks 



1 de JongJfeOOlh . The clusters' magnitude were converted 
3.29 (|Blanton fc Roweisl I2OO7I ) and multiplied by (M/ % = -9 



i n sol ar units, the lower end of the range quoted for single-metallicity populations by iMaraston 
(200 51). We also considered cluster mas ses derived from fitting of optical surface brightness profiles 
by iMcLaughlin &: van der Marell (|2005l ): these masses are p rimarily derived fr om y-band observa- 



tions (corrected for extinction) and assume (M/L)y = 2 (IMcLaug 



ilin 



cluste rs (particularly core-collapse candidates) are missing from the 
(|2005l ) database. 



2000 ). Some Milky Way 



McLaughlin &: van der Mare] 



We combined the two methods to select a final list of massive Local Group globular clusters. 
Excluding the anomalous cluster to Cen, only a handful of Milky Way clusters have masses > 
10 6 M : NGC 104, NGC 6388, NGC 6441, and NGC 6715. A further 8 Milky Way clusters 
(NGC 2808, NGC 7089, NGC 2419, NGC 6266, NGC 6273, NGC 6356, NGC 6402, NGC 6440) 
and two Magellanic Cloud clusters (NGC 1835, NGC 1916) had masses likely to be above 5x 10 5 M . 
But nearly all of the massive Local Group globular clusters belong to M31, a total of 58 objects. 
After examining the available Spitzer/IRAC data for Milky Way and LMC clusters, we eventually 
decided against using data on those galaxies' GCs in our analysis: the Milky Way cluster data had 
technical issues (long frametimes and in-place repeats) as well as only covering a single 5' x 5' field 
of view in all four IRAC bands. The LMC clusters NGC 1916 is located near a region of diffuse 
emission, making background subtraction for this clus ter extremely difficu lt. IRAC photometry of 
a larger sample of LMC clusters has been reported by lPessev et al.l (120101 ) . 



2.3. IRAC Observations and Photometry 



The M31 bulge and some of the brightest M31 clusters had individual deep Spitzer/IRAC 
observations as part of PID 3400 (PI M. Rich). These images also include a number of additional 
clusters fortuitously located near the target clusters. The data were taken with 30-second frame 
times, 22 dither positions with 4 frames per position, and so have overall depths of about 40 
minutes per sky positio n. Saturation limits for IRA C images with these frametimes are in the 
range 8 < my e g a , < 11 (jSpitzer Science Centerl l2010i ). The brightest GCs in M31 have K s ~ 11 
(IGalleti et al.ll2004l ). so saturation of the clusters was not expected to be an issue. We retrieved the 
data from the Spitzer Heritage Archive and used the post-BCD mosaics as produced by versions 
18.7 and 18.18 of the pipeline. 



-6- 





* , . # ■ / 4 

♦ ' * 

» - ;'**'•*'•-*•?' 


• ■ tx* ■ •■' 

: » 

. i ' * «■: i . 
.« v.* * • 


B225 

20 arcsec 


B178 *', •* 

L • < 
$ ■ * 


B338 

• 

••* 


B012 

# * 



Fig. 1. — IRAC 3.6/um (top) and 8.0/um (bottom) images of M31 globular clusters. All images 
are 90" square with north up and east to the left. The clusters range from K s ~ 11.5 (B225) to 
K s ~ 13.3 (B178). The left two panels show clusters in the shallow M31 mosaic (PID 3126) while 
the right two panels show clusters observed in deeper images (PID 3400). 
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The main observational dataset for the M31 clust e rs co mprises the IRAC mosaics taken as 
part of program PID 3126, described in iBarmbv et al.l (|2006l ). In brief, the observations cover a 
total area of about 3?7 x 1?0 along the major axis of M31, with an extension to include NGC 205. 
The central portion of the disk was covered to a depth of six 12-second frames per position, and the 
outer area by four 30-second fr ames per posit i on. T he mosaics used for th is analysis ar e impr oved 
versions of the mosaics used by lBarmby et al.1 (120061 ) and are described by lMould et al.l (|2008l ). Of 
the 58 clusters in our list of massive M31 globulars, 11 were imaged in all four bands in the deep 
PID3400 data, 41 were imaged only on the shallow M31 mosaics, two objects (B405, MGC-1) had 
no IRAC images available, and 3 objects (B240, B381, B383) were at the edge of the mosaic image 
and had only noisy data in 4.5 and 8.0 /iin available. Cluster B095 is located about 7" from a much 
brighter star and is contaminated by detector artifacts due to the star; we judged that obtaining 
reliable photometry in the IRAC bands was not possible. The final sample of M31 clusters analyzed 
here therefore has a total of 52 members. Sample IRAC images of M31 globular clusters are shown 
in Figure [TJ 



The typical ha lf-light radius of an M31 GC is Rh ~ 3 pc, or about 0.75" (jPeacock et al.1 120101 : 
Barmbv et al.ll2007l ). While this is smaller than the full- width at half-maximum of the IRAC point 
spread function (PSF), it is large enough that M31 GCs appear slightly resolved (see Figure [1]). 
Aperture photometry is thus more appropriate for these objects than PSF-fitting photometry. But 
how do we best account for the different PSF in the different IRAC bands, and decide which 
aperture should be used? From the point of view of capturing all of the light from a cluster, a 
large aperture is desirable. Use of 'total' magnitudes also facilitates comparison between widely- 
separated wavelength regimes, where not only the PSF but the apparent cluster structure, may 
differ (see ICohenl 120061 . for a discussion of the subtleties in aperture-matching). However, the use 
of a large aperture would result in low signal-to-noise in 5.8 and 8.0 /jm bands, where the clusters 
are faint and the background bright. Using smaller apertures has the problem that the different 
PSF results in aperture corrections that are not the same in all bands. 

To deal with the aperture issues, we took two approaches to photometry of the M31 clusters. 
To compute IRAC-only colors, we computed aperture magnitudes in a relatively small apertures, 
and dealt with the variable PSF by convolving the images in the 3.6, 4.5, and 5.8 /im bands 
to the 8.0 jim resolution. (A similar approach was used by IClemens et al.l 12009 . in their study 
of near-to-mid-infrared colors of Coma cluster galaxies.) For colors which combine the 3.6 and 
4.5 fim bands with visible-light or near-infrared bands, we measured 'total' IRAC magnitudes 
in a standard large aperture. C onvolution was done using the IRAF fconvolve task and the 
convolution kernels described by I Gordon et al.l (|2008l ) . For the clusters observed as part of the 
large mosaics, we extracted 2.5' x 2.5' 'postage stamp' images and convolved these, rather than 
attempting to convolve the original mosaics. We compared aperture photometry on the original and 
convolved deep images to show that the convolution performed as expected: the median difference 
in aperture magnitudes between convolved and unconvolved images was close to the difference 
between aperture corrections between the convolved band and the 8.0 /im. band as tabulated in 
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Fig. 2. — Distribution of Rgo — the radius containing 90% of the integrated light, in projection — 
for M31 globular clusters. 
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the IRAC Instrument Handbook (jSpitzer Science Center! l2010l ). To choose a small aperture size, 
we used the resul ts of surface brightn ess profile fitting from Hubble Space Telescope images of M31 
globular clusters (jBarmbv et al.ll2007l ) to compute Rgo, the projected radius containing 90% of the 



integrated light, for eac h cluster. Figu re [2] shows the distribution of Rgo for the 34 clusters from 



our sample analyzed by iBarmbv et al.l (|2007l ). The median of the RgQ distribution is 2'.'5, with a 
maximum of 3'.'7. These values refer to the intrinsic profile before convolution with the HST PSF; 
the effects of the Spitzer/IRAC PSF are more significant and mean that -R90 as observed in IRAC 
images will be larger. To maximize signal-to-noise and still include a significant fraction of cluster 
light, we chose one of the standard IRAC aperture sizes, a radius of 3'.'66. 

Aperture photometry was performed using the IRAF/APPHOT package. Flux was measured 
in circular apertures of radii 3'/66 (and 12 / .'2 for the 3.6 and 4.5 fim bands only), with background 
measured in an annulus of inner and outer radii 3'.' 6-8'/ 5 (14 / .'6-24'/4) subtracted. This background 
annulus for the small aperture is the smaller of two tabulated in the Instrument Handbook. We 
used it because a few clusters had nearby diffuse emission which could bias the background estimate 
in a larger annulus. We verified that in most cases the choice of background annulus had no effect 
on the resulting photometry. Calibration to the Vega system was made using the Vega IRAC-band 
flux densities given by the Spitzer Science Center. Uncertainties on the photometric measurements 
were estimated using the formula in apphot/phot, with a correction to the background noise 
term which accounts for the correlated noise introduced by pixel re-sampling during mosaicing. 
As the background variation is by far the dominant source of noise in these measurements, this 
is an important correction. Table [T] gives the photometry. We checked the internal consistency of 
the photometry by comparing measurements of the same clusters on the deep and shallow IRAC 
images. There is some scatter between the two image types, but magnitude offsets are consistent 
with zero. In the 4 IRAC bands, the offsets are (mean ± standard deviation): —0.03 ±0.04, — 0.06± 
0.08, -0.06 ± 0.07, -0.06 ± 0.09. 



2.4. 2MASS Photometry and Extinction Correction 



While globular clusters have not been extensively observed in the mid-infrared, there is a 
long history of using them to test population synthesis models in the near-infrared, so it is of 



that of 2MASS ( 


Skrutskie et al. 


2006) 


Galleti et al. 


(2004) and 


Nantais et al. 



2006), and th e 2MASS data on M 31 GCs has been examined by 



M3I GCs in the K s band based on deeper UKIRT imaging. We chose t o re- measure the M3 1 GC 
magnitudes on the 2MASS '6x' Atlas images because the observations of 
covered 60% of our cluster sample, the compila tion of[' 



Galleti et al 



Peacock et al.1 (| 20 1(1 ) only 



(2004) used the original shallow 



2MASS images of M31, and the compilation of lNantais et al.l ([20061 ) used the 2MASS Point Source 
Catalog, which might be inappropriate for the larger M31 clusters. To maintain consistency with 
the colors measured using the IRAC data, we retrieved 2MASS images of the M31 GCs and used 
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these to measure aperture magnitudes with the same large photometric and background apertures 
used for the IRAC data. We chose to use the large apertures rather than trying to match the spatial 
resolution between IRAC and 2MASS because the temporal variation in the 2MASS PSF makes 
any such matching significantly more complicated. The photometric uncertainties were estimated 
using the same procedure as for the IRAC magnitudes, which is based on the procedure detailed in 
section 8av of the 2MASS All-Sky Data Release Explanatory Supplement. TableCQgives the 2MASS 
K s photometry. We compared our aper ture photometry with the total magnitudes measured on 
UKIRT images by iPeacock et al.l (|2010l ). The UKIRT magnitudes are slightly fainter, with the 
median offset 0.08 ± 0.09 mag. This is no t too surprising, sin ce we used an aperture slightly 
larger than most of the apertures used by IPeacock et al.l (|2010l ). It could result in our 'total' 
K s magnitudes being slightly too bright; however for the purposes of computing 2MASS-to-IRAC 
colors this effect should cancel out. 

While extinction has a much smaller effect on photometry in the mid-infrared than in visible- 
light or even the near-infrared, the effect is not completely negligible. The e xact shape of the 



extinction law in the mid-infrared appears to vary wi t h position in the Galaxy (llndebetouw et al 



2005 



Roman-Zuhiga et all 120071 : iFlahertv et al.1 120071 : iNishivama et al.ll2009l : iGao et al 



200 



J but 



all of these analyses agree that the extinction at IRAC wavelengths is greater than would be 
expected from a power-law extrapolation of the near-infrared extinction curve. In t he analysis that 



follow s, we correct the cluster colors for extinction using the E(B—V) values given bv lCaldwell et al 



(|201ll ). These values include both foreground extinction a nd extinction inter n al to M31. For this 
an alysis we u s e the bandpass-averaged extinction values of Indebetouw et al. (j2005l ). as tabulated 
in lOao et all (l2009h : A 3S /A Ks = 0.56 and A\/A K , = 0.43 for the other three IRAC bands. We 
assume A K JE{B - V) = 0.35 dCambresv et alJbood ) and the standard A V /E(B - V) = 3.1. 



2.5. Population Synthesis Models 



While there a re numerous population s ynthe sis models curren tly available — two of the most 
widely-used being iBruzual k, Charlotl (|2003l ) and iMarastonl (|2005l ) — to our knowledge only three 
models give explicit predict ions for magnitudes in the IRAC bandp asses. These are GALEV 
(|Kotulla et al.ll2009l ). FSPS dConrov et a lihood: Iconrov fc Gunnlhoioh. and t he SSP models pro- 
vided by the Padova group through their STEL webservej^ ( Marigo et al. 2008 ; Girardi et al. 2010 , 
hereafter referred to as 'Padova'). While many other models can be used to generate predicted 
IRAC magnitudes from spectral energy distributions, using three model families suffices for this 
first analysis. These models contain some of the same input ingredients but differ in other ways, 
such as the treatment of thermally-pulsing asymptotic giant branch and horizontal branch stars. 
GALEV also provides chemically consistent modeling of galaxies while the other two do not, al- 
though this difference is not relevant for the present study. 



^http: //stev. oapd. inaf . it/ 
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Population synthesis models have mostly focused on the visible-light regime, so it is reasonable 
to ask what the input ingredients of these models are in the infrared. The treatment of red giant 
and asymptotic giant branch stars is expected to be particularly important, since t hese evolution- 
ary stages dom inate the light of old populations, especially in the near- infrared (jBuzzonil Il989l ; 
Marastonl 120051 . and many others). All three models were consider here use Padova isochrones, so 
an important factor in any different outputs is t he stellar libraries. GALEV and FSPS both use the 



BaSe L stellar libraries: version 2.2 for GALEV (jLejeune et al.lll998l ) and versio n 3.1 (I Wester a et al.. 
2002 ) for FSPS. BaSeL is based mainly on the stellar atmosphere models of Kurucz ( 19921 ). but 



those models do n ot inc lude the coolest stars, so BaSe L also includes model M giant atmospheres 
from lFluks et alJ (Il994h and bessell et al.l (|l989l . Il99lh. The Padova SSP mo dels use some of the 



same input ingredients — newer ATLAS models from lCastelli Kuruczl ( 20031) and thelFluks et al 
(|1994l ) models — but without going through th e BaSeL library. Since the iFluks et al.1 ()1994l ) models 



extend to wavelengths up to 12.5/mi and the lBessell et al.l (|1989l . Il99ll ) models only to 4.1/iin, we 
can infer that the former models provide the basis for the predictions tested here; however, this 
is not explicitly stated in any of the model descriptions. The thermally-pulsing asympt otic giant 
bran ch can also be an important contributor to the integrated light of star clusters (jMaraston 
20051 ). although it is not exp ected to be importan t at glo bular cluster ages. For these stars, FSPS 



uses empiric al spectra f r om Lancon &: Mouhcinel ( 2002h while the Pado va models use syn thetic 
spectra from lLoidl et al.l (|200ll ). Both the lLancon &: Mouhcinel (|2002l ) and lLoidl et al.1 (|200ll ) spec- 
tra extend only to 2.5 /im. The Padova models extend the AGB spectra beyon d 2.5 / /m with a 
Rayleig h-Jeans model: FSPS extrapolates carbon star spectra with lAringer et al.l (|2009l ) synthetic 
spectra. iKotulla et al.l (|2009) did not describe any special treatment of TP- AGB stars in GALEV. 



For al l three models, we generated predicted SSP colors in the 2MASS and IRAC bands for 
Salpeter ( 1955 ) (0.1 — 100 M ) initial mass functional and ages from 1 Gyr up to the oldest 



available. We retrieved SSP predictions from the online GALEV interface using all available fixed 
metallicity values [M/H] = log(Z/Z Q ) = (-1.7, -0.7, -0.3, 0.0, +0.3). For FSPS, we used version 
2.3 of the code to generate predicted magnitudes for simple stellar populations with metallicity 
values [M/H] = (—1.68, —0.69, —0.30, 0.0, +0.20). The default options for other parameters were 
used. For the Padova models, we used metallicities o f [M/H] = (-1.68, - 0.70, -0.3 0, 0.0, +0.20), 



with t he default parameters: evolutionary tracks from lMarigo et al.l (|2008l ). with the lGirardi et al 



(|2010l ) case A correction, and no circumstellar dust (although see below). For all models, no intrinsic 
dust extinction was applied. We interpolated the GALEV magnitudes to [M/H] = +0.2 so that all 
three models represented the same five metallicity values. These metallicity values adequately cover 
the range of metallicities for c lusters in the present sample (—1.9 < [M/H] < —0.1); although M31 
has lower- metallicity clusters (jCaldwell et al.ll201ll ). they are not present among the most massive 
objects. 



3 Comparison of all three models' predictions for ISalpeterl (1955) and iKroupal ( 20011 ) initial mass function showed 
that this choice of IMF has a negligible effect on predicted colors of old stellar populations. 
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Fig. 3. — Predicted IRAC colors from three pop ulation synthesis mo dels, as a function of sim- 



ple stellar populatio n age. Left column: F SPS IConrov et al 



((20091), center column: GALEV 



(Kot ulla et al J 12009] ), right column: Padova (IMarigo et al.ll2008l ). The same set of metallicities is 
plotted for all 3 models. 
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Figure [3] shows predicted colors from the models as functions of age. At old ages, all models 
show a rather small range of predicted mid-infrared colors, only a few tenths of a magnitude. Al- 
though the colors mostly become slightly bluer with age, the size of this effect is very small. There 
are some systematic differences in the models' predictions, with, in general, the FSPS models pre- 
dicting the greatest color variations over the displayed age range. While there is general agreement 
that higher metallicities result in redder colors within a model family, the different models do not 
agree on either the size of this effect or its location. For example, the FSPS K s — [3.6] color for 
the most metal-poor model is redder than the same color for the most metal-rich Padova model. 
Each of the models has at least one metallicity track that behaves differently from others in the 
same family, but which metallicity is the outlier changes depending on the color and model family. 
Since all three models use the same set of isochrones, the difference in predicted colors must be 
due to some other factor. One likely contributor is the different spectral libraries used for post- 
main-sequence stars. Another could be numerical issues surrounding the s tability of integration 



over short-lived, but luminous stages (jMaraston 



2005 



Tinslev k GumJll97a ). 



While Milky Way glo bular clusters have not shown evidence for intracluster dust (jMatsunaga et al 



2008 



Barmbv et al.ll2009l ). mid-infrared obser yations have provided important evide nce for dust in 



the envelopes of individual mass-losing stars (jBover et al 



2010 



Origlia et al.ll2010l ). So it is rea 



sonable to ask what the effects of dusty mass loss would be on the integrated mid-infrared colors 
of a clusters. The Padova models can optionally include these effects for stars wit h significant 
mass loss: red supergiants, TP-AGB, and upper-RGB stars (see iMarigo et al.l 120081 . for details). 
Figure [4] shows the result of adding circumstellar dust to the Padova simple stellar populations. 
This figure shows the Padova models with three different prescriptions for cir cumstellar dust: sil- 
icate dust for M stars and carbon-rich dust for C stars (IBressan et al.lll998l ). 100% A10 x for M 
stars and 100% amorphous carbon for C stars, and 60% silicate /40% AlCy for M s tars and 85% 
amorphous carbon/15% SiC for C stars (the latter two both from lGroenewegenlhood ). Comparing 
with the right panel of Figure [3] shows that including circumstellar dust emission causes the colors 
to vary more with age (all of the Padova models are plotted with the same time resolution) and 
become redder, particularly at the longer wavelengths. These effects are stronger for ages < 8 Gyr: 
models at the oldest ages are close in color to the non-circumstellar-dust version. In most cases, 
the choice of circumstellar dust prescrip tion does not have a large effect on the predicted colors. 
The exception is [5.8] — [8.0] where the IBressan et al.l (|1998l ) prescriptio n results in b l uer co lors, 
particularly at lower metallicity, than either of the two prescriptions from iGroenewegenl (|2006l ) . To 
maintain compatibility with the other two families, the dusty circumstellar models are not shown 
in the comparison with the cluster data, but these effects will be considered in the discussion. 
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3. Analysis and Discussion 

3.1. Cluster color distributions 

Before comparing to models, the most basic analysis of the globular cluster colors is to examine 
their distribution. Figure [5] shows the distribution of the four colors constructed with bandpasses 
adjacent in wavelength. Overplotted are Gaussians normalized to the same area, with the data- 
derived means and standard deviations: K s — [3.6] = 0.13 ± 0.11, [3.6] — [4.5] = — 0.04 ± 0.05, 



[4.51 - [5.8] = 0.09 ± 0.07 and [5.8] - [8.0] = 0.09 ± 0.12 (mean ± standard deviation). IPahre et al 



(J200J) gave theoretical colors of M0 III stars in the IRAC bands of [3.6] -[4.5] = -0.15, [4.5] -[5.8] = 
+0.11 and [5.8] — [8.0] = +0.04. The mean colors of the cluster sample are formally consistent with 
the naive expectation of zero magnitudes, and in the longer wavelengths with the theoretical colors. 
The mean [3.6] — [4.5] color of the clusters is slightly redder than the expected giant color, as are 
the model colo rs described in Se ction 13.21 so we suggest that there may be a problem with the 



value given by IPahre et al .1 (|2004l ). The [3.6] — [4.5] color shows the narrowest color distribution, 
unsurprising since these two bands have the most similar PSF and do not suffer from any possible 
aperture-matching problems between 2MASS and IRAC data. At the longer wavelengths, where 
the clusters are fainter and the background brighter, the observational uncertainties are larger. 
Contributions by diffuse PAH emission could also skew the 5.8 and 8.0 (im fluxes, although this 
should be minimized by the close-in background aperture used for photometry. 

As Figure[3]shows, the expected range of mid-infrared colors for old stellar populations is rather 
small, a few tenths of a magnitude. Are the distributions that we measure the result of a single 
intrinsic color, broadened by observational uncertainty? One argument against this is that, for the 
IRAC-only colors, the standard deviations of the color distributions are broader than the average 
photometric uncertainties, by factors of 1.5-2.5. Applying the DAgostino- Pearson normality test 
to all four color distributions, we can reject the hypothesis of normality in the [3.6] — [4.5] and 
[5.8] — [8.0] colors, but not in K s — [3.6] or [4.5] — [5.8]. We also compared the color distributions 
to Gaussians with zero mean and standard deviation given by the average photometric uncertainty 
using a Kolmorgorov-Smirnov test. In all cases the hypothesis that the observed distributions were 
drawn from the zero- mean distribution was rejected. We conclude that the cluster colors have some 
intrinsic spread and that the observed distributions are not merely delta functions broadened by 
observational uncertainty. Mis-estimates of reddening should have only a small effect on colors 
derived from adjacent band and are unlikely to cause significant scatter. A spread in ages is a 
possible reason for a color distribution spread; this is discussed further in Section 13.21 Another 
possible explanation is that the relatively poor spatial resolution of IRAC leads to contamination 
of some clusters' colors by sources nearby in the M31 disk. 




Fig. 5. — Distribution of near- and mid- infrared extinction-corrected colors for massive M31 globu- 
lar clusters. Gaussians with mean and standard deviation computed from the data are overplotted 
for reference. 
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Fig. 6. — Near and mid-infrared colors (corrected for extinction) as a function of metallicity for 
massive M31 globular clusters (points) and SSP model predictions (lines). Each of the model 
prediction lines represents one model family at a specific age between 2 and 12.5 Gyr. Solid 
black symbols are individual M31 clusters; open symbols (red in the online version) represent the 
luminosity-weighted mean colors of clusters in quintiles of metallicity. 
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3.2. Observed and model colors 



While, by design, the M31 clusters in our sample cover a fairly small range in mass, they do 
cover a range of about 2 dex in metallicity. Comparing models to cluster colors as a function of 
metallicity allows us to see how the ranges of the data and models coincide. Figure [6] shows these 
comparisons for the same colors as shown in Figure [5j The models plotted are a subset of those 
shown in Figure El chosen to represent the range of model colors for ages of 2, 8, an d 12.5 Gyr. For 



Caldwell et al 



metal licity values of the clusters, we used the spectroscopic [Fe/H] measurem ents of 
(hoilh and converte d to total metallicity with [M/H] = [Fe/H] + 0.94[a/Fe] (|Thomas et al.ll2003h 



Beasley et al.l (120051 ) found a large spread in the degree of a-element enhancement for M31 clusters; 
we use a typical value of [a/Fe] = 0.2. As well as data for individual clusters, the Figure also shows 
the luminosity-weighted mean colors of clusters grouped in quintiles of metallicity. In general, the 
clusters cover a wider range of color than do the models, but the models typically fall within the 
range of the cluster data. The reddest model in the K s — [3.6] and [3.6] — [4.5] plots is the FSPS 
2 Gyr model, which does not appear to be a good fit to the data; however, the 2 Gyr Padova model 
does match the data in K s — [3.6] , as does the 2 Gyr GALEV model in [3.6] — [4.5] . The observational 
scatter in these two colors is large enough that older models are also broadly consistent with the 
data. 

At the longer IRAC wavelengths, the range of model colors is very small and the models are 
generally in good agreement with each other. In [4.5] — [5.8], the cluster mean colors are slightly 
redder than all of the models by up to 0.1 mag, but the offset is statistically significant in only 
two of the five metallicity bins. The scatter of some objects to redder colors could be due to 
systematic effects in the photometry (incomplete background subtraction, failure to completely 
match apertures with the convolution method), or, more interestingly, it could indicate that some 
clusters' colors are affected by circumstellar dust emission. In [5.8] — [8.0], the mean cluster colors 
are generally in agreement with the models except possibly at the highest metallicity. For the most 
part, the redder colors predicted by the Padova models for stars with circumstellar dust do not 
appear to match the cluster data. This lends support to the idea that these clusters are truly old 
objects and their infrared light is not dominated by emission from circumstellar dust. 

Optical to near-infrared colors (e.g., V — K s ) are good metallicity indicators in old stellar 
populations, because they are sensitive to the temperature of the red giant branch. Optical-to- 
mid-infrared colors have been less well studied. ISpitler et al.1 (J2008J) measured V — [3.6] for GCs in 
NGC 5128 and NGC 4594 and found that this color showed a good correlation with spectroscopic 
metallicity. To see if the same is true for M31 GCs, we compute d optical-to-mid-infrared colors 
using the F-band photometry t abulated by ICaldwell et al.l (|2009i ). corrected for extinction using 
the E(B — V) values given by ICaldwell et al.l (|201ll ). The se colors are much m ore sensitive to 
reddening than the infrared-only colors considered above. ICaldwell et al.l (|201ll ) estimated the 
typical reddening uncertainty to be 0.1 mag: changing a reddening estimate by this much would 
correspond to a change of nearly 0.3 mag in reddening-corrected V — K s . 



-19- 




Fig. 7. — Optical to near /mid-infrared colors (corrected for extinction) as a function of metallicity 
for massive M31 globular clusters (points) and SSP model predictions (thick lines). Thin lines in 



Cohen et al 



the to p two panels are fits of V — K s against metallicity for Milky Way clusters from 
(2007) and V - [3.6] against metallicity for GCs in NGC 5128 and NGC 4594 from ISpitler et al 
(|2008l ). Symbols as in Figure error bars on the data points do not include uncertainties due 
to reddening correction, which could be up to 0.3 mag. The very blue cluster at JM/H] = —0.6 
is B037, a cluster with a large (E(B — V) > 1.3) and highly uncertain reddening (ICaldwell et al 



2011 



Barmbv et al J 12002 ). 
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Figure [7] shows the V — K s , V — [3.6] and V — [4.5] colors as a function of metallicity for M31 
GCs. Overplotted are the relations predicted by the GALEV, FSPS, and Padova models (ages 2 
and 12 Gyr), as well as fits to the color-metallic ity relationships for other GC sys t ems: V — K s 
for Milky Way clusters from bohen et al.l (|2003) and V - [3.6] from ISpitler et al.l (feood ). In all 
three model families, older models are redder at a given metallicity than younger models. There 
is a significant spread in predicted color with metallicity: the GALEV models are almost always 
redder than the Padova models at a given metallicity, with the FSPS models intermediate. The 
Figure shows that the models do not predict a strictly linear color-metallicity relation for any of 
the colors; particularly for V — [3.6] and V — [4.5], all of the models show a slight change in slope at 
[M/H] = —0 .7. The empirica l fit in V — K s captures the model behavior reasonably well, although 
as noted by IConrov fe Gunnl (|2010i ). the models are redder than the data at high metallicities. 
In V — [3.6], the fit matches the model s well at high metallicities but is bluer than the models 
at low metallicities. ISpitler et al.l (|2008l ) found the same result, although interestingly they were 
comparing to a completely different set of models. 

The cluster data points plotted in Figure [7] are bluer than the models except at the lowest 
metallicities. A problem with the photometry does not provide an immediate explanation. Since 
our large-aperture IR magnitudes could slightly overestimate the clusters' brightness (because of 
contamination by nearby objects), our photometry seems more prone to make the clusters too red, 
rather than too blue. A number of other systematic effects could come into play: 



One possibility is that our method of using the average a-element abundances to compute 
[M/H] is too simplistic, although the [M/H] values would have to change drastically to make 
the data points compatible with the models. 

Another possibility is that the clusters are much younger than expected and should therefore 
be compared to younger models. We think this is unlikely: deriving ages from mid-IR model 
predictions and colors alone (clearly a risky idea given the limited age sensitivity shown in 
Figure [3]) would suggest ages younger than 2 Gyr. Such young ages would have a substantial 
effect on optical and near-infrared colors, and there are no suggestions from recent spec 



troscopic or SED -fitting studies of M31 clusters (e.g. iRev et al. 



2007; Caldwell et al. 2009 



Wang et al.ll2010l ) that the brightest clusters are all so young. 



3. Reddening corrections could play a role: the metal-rich clusters have a higher avera g e red - 
dening than the metal-poor clusters, so if the E(B — V) estimates of I Caldwell et al.l (|201ll ) 
were too high, the corrected colors of metal-rich clusters would be more strongly affected. As 
an example, we found that reducing the reddening values of all clusters by 20% brought the 
data points and models into better agreement over the entire metallicity range@ 



4 The reddening law itself is unlikely to be the problem, since almost all of the correction to the V— infrared colors 
is in Ay and there is no strong evidence that M31 has a different value of Rv- 
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4. All three model families could have some common flaw in the spectral libraries that causes 
them to predict too much infrared flux, or too little U-band. Some SSP models are known 
to have problems reproducing globular cluster colors in the infrared (jConrov &: Gunnll2010l . 
Figure 6), so a problem with the less- well-calibrated IRAC bands seems plausible. 



5. Finally, the results of lStrader et al.1 (|2009l ) on the mass-to- light ratios for M31 clusters suggest 
an interesting alternative. Those authors found M/L to be lower for metal-rich M31 clusters 
than for metal-poor clusters, contrary to predictions from SSP models. While younger pop- 
ulations are expected to have lower M/L, IStrader et al.1 argued against the metal-rich M31 
clusters being even intermediate-age, and concluded that variations in the cluster mass func- 
tion, either initially or through dynamical evolution, were more likely explanations. A flatter 
mass function, required to match the mass-to-light ratios, should qualitatively also result in 
bluer predicted colors: fewer low-mass stars would mean less U-band light (dominated by 
the main sequence), and more high-mass stars should mean more infrared light (dominated 
by giants). The different initial mass functions available in the SSP models do not result in 
large differences in predicted colors, but perhaps more drastic alterations to the mass function 
would change colors more dramatically. 



We consider the last three of these options to be more plausible explanations of the offset. Distin- 
guishing between them requires better-quality data than are currently available. 

If we assume that the model-data offset in V— infrared colors is due to a systematic problem 
in the models, we can still use our data to explore the metallicity sensitivity of these colors. We 
performed least-squares fits of color as a function of metallicity, weighted by the color uncertainty. 
For the clusters binned in metallicity, we found slopes i n the range 0.42-0.46 for all three colors, 
very close to the value of 0.409 derived by lCohen et al.1 (120071 ) for V — K s in Milky Way globular 
clusters. For the full, unbinned M31 sample we derive a slope of 0.39 for V — K s . For this sample 
we derive somewhat steeper slopes of 0.68 for V — [3.6] and 0.59 for V — [4.5], still shallower than 
the slope of 0.90 derived by ISpitler et al.l (|2008l ). although those authors were fitting the reverse 
function (metallicity as a function of color). The steeper slopes for the unbinned data agree better 
with the general trend of the models. Of course, there i s no guarantee th at a linear fit is the correct 
representation of the actual physical situation (and see lYoon et al.ll2006l ). but the scatter visible in 
Figure [7] indicates that a more complicated fit is not likely to be meaningful. 

The slopes of the linear fits indicate that V — IRAC colors could be a useful metallicity indicator 
for simple stellar populations. The low infrared background seen by Spitzer means that IRAC 
observations to comparable depth can be made much more quickly than ground-based near-infrared 
observations; however compared to V — K s , V — IRAC colors have the disadvantage of being 
somewh at more sensit i ve to r eddening uncertainties and having (typically) poorer spatial resolution. 
Echoing ISpitler et al.l (J2008J), we caution that these colors are not well-calibrated for super-solar 
or very low ([M/H] < —1.9) metallicities. Using the M31 clusters to calibrate the relation would 
require careful matching of apertures in photometric bands differing by almost a factor of 10 in 
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wavelength (see ICohen et al.l 120071 . for a description of the difficulties involved in matching just 
between visible and near-infrared) and is not a trivial task even for these bright clusters. Using 
more distant, unresolved clusters to calibrate the visible-to-IRAC color /metallicity relation might 
be a better approach. 



4. Summary 

We have presented mid-infrared photometry of massive globular clusters in M31, in order to 
test the color predictions of population synthesis models. To our knowledge, this is the first such 
test of the models in all four IRAC bands. The three different model families discussed here differ 
somewhat in their predicted colors and in the systematic effects of age and metallicity. As might be 
expected, models including circumstellar dust emission predict redder colors, particularly at longer 
wavelengths, and more color variability with age, particularly for ages < 8 Gyr. Except for the 
[5.8] — [8.0] color, the exact choice of circumstellar dust prescription does not make a substantial 
difference. The present level of observational precision does not permit distinguishing between 
models; the considerable observational uncertainties, especially at the longer wavelengths, are due 
to the bright background from the galaxy and the falling cluster spectral energy distribution. 
Within these uncertainties, the data and models are mostly in agreement. We find a small but not 
statistically significant offset between models and data in the [4.5] — [5.8] color; if real, this could 
be an indication of emission from circumstellar dust in some clusters. Finally, we find that visible- 
to-IRAC colors are bluer than model predictions, for reasons which are unclear. These colors could 
still be feasible metallicity indicators for simple stellar populations but further work is required to 
establish accurate calibrations. 

Facilities: Spitzer (IRAC), FLWO:2MASS () 
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Table 1. IRAC photometry for massive M31 globular clusters 
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.01 


11 


.53 


± 


0. 


01 


-1 


1 





.15 


B358 


13 


33 


± 


0. 


.01 


13 


.32 


± 





.01 


13 


21 


± 





01 


13 


.19 


± 


0. 


02 


12 


.96 


± 





12 


12 


.76 


± 





03 


12 


.74 


zt 


0. 


02 


-1 


.9 





.06 


B373 


12 


.68 


± 


0. 


.01 


12 


72 


± 





.01 


12 


.71 


± 





02 


12 


.71 


± 


0. 


07 


12 


20 


± 





.03 


12 


.18 


zt 





01 


12 


2(1 


± 


0. 


01 


-0 


.5 





.25 


B386 


13 


.15 


± 


0. 


.01 


13 


.17 


± 





.01 


13 


.12 


± 





02 


13 


00 


± 


0. 


03 


12 


.72 


± 





11 


12 


.44 


zt 





01 


12 


.43 


zt 


0. 


01 


-1 


1 





.18 


B472 


12 


.85 


± 


0. 


(12 


12 


.80 


± 





.02 


12 


.73 


± 





.03 


12 


.46 


± 


0. 


04 


12 


39 


± 





.10 


12 


.20 


zt 





03 


12 


13 


± 


0. 


03 


-1 


.0 





12 


G001 


11 


.34 


± 


0. 


.01 


11 


.37 


± 





.01 


11 


.22 


± 





01 


11 


.12 


± 


0. 


01 


10 


.61 


± 





.02 


10 


.55 


zt 





01 


10 


57 


± 


0. 


01 


-1 


.0 





.10 


MITA140 


12 


.72 


± 


0. 


02 


12 


.75 


± 





.01 


12 


.68 


± 


0. 


02 


12 


.56 


± 


0. 


05 


12 


40 


± 





13 


12 


30 


zt 





05 


12 


.23 


± 


0. 


04 


-0 


2 





.86 



Note. — Values are not corrected for extinction. Measurements in the 5.8 and 8.0 fim bands arc made only in the smaller aperture because of the 
brighter background at these wavelengths. 

a Naming convention is that of Gallcti ct al 

b From [Cald well et ah! j201l|) . 

c This cluster is located near a much brighter star; adequate large-aperture photometry was not possible. 
These clusters are located near the bright bulge of M31; adequate large-aperture photometry was not possible. 



